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0'\cX/waygiA;  0  ,  abstract 

A   recurring   problem   in   solid   state   physics   is   that   of 
determining   the    electrostatic   field   or   potential  at   a   point 
within  the   crystal   lattice   arising  from   the   molecules   of   the 
lattice  themselves.      A  typical   problem   is    that   of   finding 
the   electrostatic   binding   potential   of   an  ionic   lattice, 
assuming   point    changes. 

Nijboer,    de  Wette   and   others   have   worked   toward 
perfecting   a   method  for    calculating   the    electric   field 
arising  from   a   lattice   of   dipoles.       This   work   culminated 
in   a   paper   by   de  Wette   and   Schacher   detailing  a   procedure 
which    ensures    a   rapidly   converging   summation  of   the   three 
dimensional   lattice   sums   involved. 

The   present   paper   is   an    extension   of   the   latter, 
presenting   a   computer    program    for    evaluating   dipole   fields 
in   a   crystal  lattice   using   the   equations   obtained  by   de  Wette 
and   Schacher. 
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BACKGROUND 
A   recurring    problem    in   solid    state    physics    is    that    of 
determining   the   electrostatic   field   or   potential   at   a   point 
within   the    crystal   lattice   arising   from    the    molecules    of   the 
lattice   themselves.      A   typical   problem   is    that   of   finding 
the   electrostatic   binding   potential   of   an  ionic   lattice, 

assuming   point   changes. 

1    2 
Nijboer,    de  Wette    '      and   others   have   worked   toward 

perfecting   a   method   for    calculating  the   electric   field 

arising   from   a   lattice   of   dipoles.       This    work   culminated 

in   a  paper   by   de  Wette   and   Schacher^    detailing  a   procedure 

which    ensures   a   rapidly   converging   summation  of   the   three 

dimensional   lattice   sums   involved. 

The   present   paper   is   an   extension   of   the   latter, 

presenting  a   computer   program   for    evaluating   dipole   fields 

in  a   crystal   lattice   using   the   equations   obtained  by   de  Wette 

and   Schacher. 


B.R.A.    Nijboer    and    F.W.    de  Wette,    Physica    14,    422 
(1958) 

2F.W.    de   Wette,    Phys.    Rev.     123,    103    (1961) 

3 
F.W.    de   Wette   and   G.E.    Schacher,    Phys.    Rev.     137, 

A78    (1965) 


THEORY2* 
Virtually   all   problems    concerning   crystal   self   fields 
reduce   ultimately   to    evaluating  a   lattice   sum   of    the 
general  form: 

S  (X.lLJD  =  L~  fCgrOOf ,  yrbOx  M-CO*) 

where   x>    y»    and    z   are   coordinates    of   the   point   at   "which 
the   sum   is   to   be    evaluated;    f    (     )    represents   the   contri- 
bution to    the   sum   from   a   particular   lattice  point;    and 
aa_]_,    bag,    cao    are   the   coordinates    of   the   lattice   points 
■whose   contribution  is   being   considered    (a-^,    a,,    and   ao 
being   lattice  basis   vectors). 

Finding   the   electric   field   due   to   a   lattice   of    dipoles 
reduces   to    evaluating   such   a   sum,    where  the   summation 
is,    mathematically,    conditionally   convergent.       That   is,    the 
value   of    the   sum   depends   on   the   order   in   which   the   sum- 
mation is   carried   out    ( analagous   to   the  fact   that   the 
internal   field   of   a   crystal   depends   on  the   shape   of   the 
crystal).       Thus,    any   procedure   for    evaluating   dipole   sums 
must,    at   some   point,    make   an  assumption  of   a   particular 
crystal   shape. 

For   a   complete   description   see    (2) 
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In   the   procedure   used  here   these   conditionally 
convergent    triple   sums   over    the   real   lattice   are   converted 
into   absolutely    (rapidly)    convergent   double   sums    over   the 
recipricol   lattice.       Converting  from   a   triple   to   a   double 
sum   forces   an   assumption   of   a   slab-shaped    crystal   with 
the   short   dimension  in   the    z-direction    (perpendicular   to 
the  plane   described  by  the   lattice   vectors   a.    and   ap  • 
See   figure    1    for   relations   between  x,    y,    z    coordinates 
and   the   lattice   basis   vectors). 

To   achieve   the   conversion  to   an  absolutely   convergent 
double   sum,    the   original   sum   is   first   Fourier   transformed, 
making   use   of   the   property   of   lattice   sums    that   the   sum- 
mation of   a   function   over   a   lattice   in  real   space   is    equal 
to    the   sum    of   its    Fourier   transform    over   the   recipricol 
lattice    (in   Fourier    space).       In   Fourier   space   it   is   then 
possible   to    carry   out   the   z-summation  analytically,    leaving 
an  absolutely   convergent   double   sum   in   recipricol   space 
(with   the   assumption   of   a   slab-shaped   form). 

A.  distinction   is    made  between   the   cases    where  the 
evaluation  point   does    or   does   not   lie   in   an   x-y   plane   of 
dipoles    (z  =  cao    or    z^ca,,    c   an  integer).       In   the  first 
case   the  function   to   be   summed   has   a   simple   pole    (with 
respect   to   the    z-summation)    at    z  =  0 .       In   this    event   an 
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auxilliary   convergence  function  is   used  and   the   sum   is 
broken  up  into   two    sums,    one  in  real   space   over   the   x-y 
plane   containing  the    evaluation   point   and   one  in  recipricol 
space   over   the   rest   of   the  lattice.      In  the   second   case 
the   function   is    everywhere   regular  and  is   Fourier   trans- 
formed  as   it   stands. 

In   either    event   the   method  is   directly  applicable   only 
to   lattices   in   which   all   dipoles   are   of   the   same   magnitude 
and   orientation.       To    evaluate   the  field   of   a   lattice   contain- 
ing  more   than   one   orientation   and/or    magnitude   of   dipoles 
it  is   necessary   to   divide   the   lattice  into    sublattices,    each 
containing   only   one   kind   of   dipole,    evaluate   the  field  for 
each   sublattice,    and   combine  the  fields   by   superposition. 

The   summation   procedure,    applied   to    a   sublattice, 
produces   nine  numbers    S-.    :    i,j   =   x,y,z    ;    where   Sm   is   the 

J  J 

i-component   of   the   electric   field    (in   Coulombs)    due  to   a 
lattice   of   unit   dipoles    (one   electronic    charge   and   length 
1A)    oriented   in  the  j-direction. 

While  program   inputs   are  in   terms   of   the   normal 
crystallographic    quantities    (a.,,    ao,    a-  i0(,)S)  7fl  $     )    the 
calculations   and   outputs   are  in   terms    of   a   cartesian 
reference   frame  related   to   the   crystal   basis   vectors   as 
detailed   in  figures    1-3  • 
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Figure   1:      Unit   Cell:       1   dipole  per   unit   cell  volume, 
Basis   vectors   a^,    a^,    a.o.      Crystal  angles:   <x*Z.fl*,.fl*    ; 
/3-Z^>£  ;   X*/-£L,CU  >  S*Z.JLIQ9 

x  axis   parallel  to   aj_,    y   axis   in  a^   -   a^>   plane  perpendicular 
to   x  axis,    z   axis   perpendicular   to   x-y  plane. 


Figure  2:      Point   at   which  field  is   to  be   evaluated: 
located  with  respect   to   lattice   point  by  fractions   of 
basis   vectors^  *j,  a, +JjlGU+JsO$  J   °^Ji^/      - 


Figure  3?      Decomposition  of  arbitrary  dipole! 


Then  given   a   lattice   of   arbitrary  dipoles  £,    as   in  figure 
3,    we  have: 

px   =   p   sin  9    cos  0 
p      =   p   sin   0  sin  0 

J 

pz  =  p  cos  a 

so   that   the  field  at   some   point  j(j^,    jg,    jo)    is   given  by: 

Ex(j>    =   PxSxx(J)+PySxy(J)+PzSxz<J) 
Ey(j)    =   PxSyx(j)+pySyy(j)+pzSyz(j) 
E2(j)    =   PxSzx(j)+pySzy(j)+p2Szz(j) 
The   nine   S-   are  not  all  independent   of    course.      They 
satisfy 

Sjj   =   Sji 

sxx  +  syy  +  szz  "  ° 
for   all   lattices.      Other   special   relations    (such  as    S..   =   0, 

i   4  j»    for    cubic   lattices)    are  also   obtained   for   lattices 

with   a  high   order   of    symmetry. 

In   a   lattice   containing    more   than  one   kind   of   dipole 

two   special   considerations   arise  in  decomposing  the   lattice 

into    sublattices.       First   and   most   obviously,    one   must 

choose   the  j's   in   each   sublattice,    determining   the   point   at 

■which  the  fields   are   to   be    evaluated,    such   that   they   refer 

to   a   common  point   in  the   composite   lattice.       The    second 

consideration  refers   back  to   the   shape   dependence   of   the 
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sum.      In   computing   the   sums   over   each   sublattice,    the 
shape   dependence   has   introduced   a   depolarization  factor 

O   Q 

L    =   i+7r/(vol   of   unit    cell   in  AJ )    into   the   S         sum.       Then 

zz 

to    superpose  the   fields   from    several   sublattices    with 
various    z-directions    (referred  to   the   composite   lattice) 
it   is   necessary   to    modify  all   the   Sjj   to   bring  all   sublattices 
to   a   common  x-y-z   orientation  by   manipulating   the 
depolarization  factors   for    each   sublattice.      A   similar 
consideration   arises   in   comparing   fields    obtained  by   this 
method   and   those   obtained  by   some   other    method   which 
assumes    a   different    crystal   shape   or   with    experimental 
results   on  a   specific    crystal   of   known   shape.       Kittel 
considers   the  general   problem   of    depolarization  factors 
in  some   detail   and  tabulates   the   depolarization  factors 
for   most   commonly   used   shapes    (spherical,    cubic,    slab, 
etc.)      An   example   of   such   considerations   for   a  benzene 
model   is    contained  in  Appendix    D. 


5 
C.    Kittel,    Introduction  to    Solid   State   Physics, 

2nd   edition,    158,     (John  Wiley  and   Sons,    Inc.,    New    York, 

1956) 


PROGRAM    COMPOSITION 

The   actual   translation   of   the   method   of   Schacher 
and   de  Wtftte  is   relatively   straightforward  although  quite 
lengthy. 

The   program   is    -written  in  the   FORTRAN   60   language 
for   the   CDC    1604   computer   and  is   compatable   with   most 
FORTRAN    monitors.       The   program   is   divided  into   a   main 
program,    containing  all   input-output   routines   and  a   segment 
to    choose  a   subroutine,    and   eight   subroutines   to   do   the 
actual   summation. 

Simplified   flow    charts   of   the   entire   program   and 
various    segments   are   shown  in  figures   2+-7. 


6 
Appendix  A   contains   the   program    listing  and   details 

of   timing,    storage,    etc. 
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Figure  k :       Flow   Main   Program   -   LATSUM 


START 


( 


I 


nput    lattice   parameters    and   convergence 
limits  


I 


3 


(Output   lattice   parameters^ 

"  i 

if    error      T7Z 7    ,    ..  . 7 ; 7 1 

___ |Test   lattice   parameters   for    correctness   1 

in   input  H 

,_..._ . Jt 

Compute  basic    constants 

l 


Examine   lattice   parameters    to   determine 
which   subroutine   to   use      


/Cutput    error\ 
y     code J 


Do   actual 
summation 


Subroutine   SUM1 


p"Hsubroutine   SUM2J 


Subroutine   SUM3 


A      Subroutines    SUM4,SUM6, 
GAM1,  GAM2 


*> 


4 


_  Subroutines;  SUM5 ,  SU M6 , 
ni   a  JniocT  eob    ■  :,'i         LOAMl ,  GAM2 


>nJBlalSA 


Output   Sjj's  --  A    ) 


(Output    subroutine   used) 

I        /      I  

(Output-^ results   of ,-:  convergence;  tests) 


END 


r 


v  v 

89  Y 


f8TAJ   C 


SI 
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Figure   5:      Common   Flow-Subroutine   SUM1,    SUM2,    SUM3 

Main  program   uses      SUM1   if   XU3^0,    XDEL^O 

SUM2  XJ3/0,    XDEL=0,XJl+XJ2/0 

SUM3  XvJ3^0,    XDEL=0,XUl+XvJ2=0 


ENTER    FROM 
L AT SUM 


Compute  basic   constants 


Compute   contribution  to   S 
for   lattice  points   on 
Al   and  A2   axis 


y 


No 


Compute   contribution   to    S 
for   lattice  points   in 


A1-A2   plane 


ij 


No 


Yes 


RETURN 

TO   LATSUM 
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Figure  6:       Common   Flow   -   Subroutine   SUM4,    SUM5,    GAM1 

Main  program   uses    SUM4  &   SUM6   if   XU3=0,    XDEL/0 

SUM5    &   SUM6    if    XU3=0,    XDEL=0 


SUM4   or   SUM5 


ENTER    FROM 

LATSUM 

1 


Compute   basic 
constants 


Compute   recip.  lattice     — 
contribution  to    S--  4 

for   lattice   points   on 

Al   &  A2   axis 


No 


Convergent 
tests    met' 


Yes 


Compute   recip.  lattice 
contribution   to    S^ 
for   lattice  points   in 
A1-A2    plane 


No 


Subroutine 
GAM1 


^Compute   auxilliary 
convergence   functions 


RETURN  TO  LATSUM 
AND  THEN  ENTER 
SUBROUTINE  SUM6 
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Figure   7:      Flow    Subroutine   SUM6,    GAM2 
Used   whenever   XvJ3=0 
SUM6 


ENTER    FROM 

LATSUM 


Compute  basic    constants  I 


No 


Add   constant   term   to   Sjj 
1 


Compute   real   lattice   contribution 
to    Sj:   for   lattice   points    on  Al 
— : fr.A?.   axis 


ISfiJ 


Yes 


Compute   real   lattice   contribution 
to    S-.   for   lattice   points   in  A1-A2 
plane 


If* 


Subroutine  GAM2 


Compute  auxilliary 
convergence  functions 

V(3/z )  rrx  z>  ynd 

rm>7rxz) 


No 


RETURN    TO 
LATSUM 
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As   can  be   seen   from   the   flow    charts,    a   deck   can  be 
made   up   for   a   particular   lattice   sum   by  appending   to    the 
main   program    either    subroutines    SUM1,    SUM2,    and   SUM3 
or   subroutines    SUM4,    SUM5,    SUM6,    GAM1   and  GAM2 
depending   on   whether    U    is   not    zero   or    zero   respectively. 
Further,    if    sums   are   to   be   computed   for    several    sets   of 
data   the    main  program    deck   can  be   modified  by  looping 
from   the  input   segment   to   the    END    card. 

Inputs   to   the   program    are   given   in   Table    1. 

Output   is    scattered   throughout   the    main  program. 
A   check   output   of   the   input   parameters   always   appears. 
If   there  is   an   error   in   the   input   parameters   the  following 
is   printed   and   control   returned  to   the   monitor: 

PROGRAM   STOPPED    BY    INPUT    ERROR 

DESIGNATED    BY    K= 
The  individual    error    codes   are: 

K  ERROR 

1  Al^O. 

2  A2  ^  0  . 
4  A3  ^  0  . 

10  XJ1  <0.    or    XJ1  >1. 

20  X<J2<0.    or   XJ2>1. 

1+0  XJ3<0.    or   XJ3>1. 

100  XALPH^O.    or    XALPH>90. 

200  XBET^O.    or   XBET>90. 

400  XGAM^O.    or    XGAM^180. 

1000  XDEL<0.    or    XDEL>90. 

2000  N^0 

15 


Table    1:      Program    Inputs 


LABEL 


UNITS 


DESCRIPTION 


ALLOWED 
LIMITS 


Al 
A2 

Al 


A 


Length   of   basis   vector   a^ 


aj>  0 


XJT1 
XJ2 
X<J3 


Dimen 

sion 

less 


jj    locating  pt   at   which  field 
is   to   be   evaluated  in   terms   of 
aj     (see  figure  2) 


3^Ji  *'l 


XALPH 
XBET 

XGAM 


Degrees 

and 
decimals 


Crystal   angles 


(see  fig   1) 


D  <XALPH^  90 
3<XBET^90 
3  <  XGAM  <  180 


XDEL 


AEX 
REX 
AEY 
REY 

AEZ 
RE2 


degrees 

and 
deci 


malsent 


A^J7 

R=fr 
dimen 
sion 
less 


Inclination  of  A3   axis   from 
perpendicular     .(see  figure    1) 
Enter   value  if   known   else 

er   XDEL=90   and   program 
will    eompnhfi 


D^  XDEL^90 


Error   limits   for    convergence 
bests.       Summation  will   contin- 
ue  until   the  next   term   to   be 
summed  is   less   than 
(x)  (x) 

AE(y)+RE(y)times   absolute 

(z)  (z) 

value   of   sum   at   that   time 
AEX&REX    control   Six&Sxi 
AEY&REY    control    Sjy^S   j 

AE2&RE2    control    S?_6-S    , 

iz       zi 

\  fuller   discussion  is   given  in 

Appendix  B       


NONE 


N 


dimen 

sion 

less 


ummation   cutoff.       Stops    sum- 
nation  after    covering   2N   by  2N 
;erms   in  x-y   plane   if    converg- 
5nce   tests    not    met   before   then 


N>  0 
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Multiple   errors    result   in   a   K   value   which   is   the   sum   of   the 
individual   K- values   involved. 

A   second   error    output    (and  transfer   back   of    control) 
occurs   if    N   is   so    small   that   it   is   impossible  to    check   the 
sums   for    convergence.       This   subject   is    covered   more 
completely   in  appendix  B.       The   printout   is    self-explanatory: 

PROGRAM   STOPPED    BECAUSE    N    TOO    SMALL. 

N    MUST    BE   AT   LEAST +    1    TO    ENSURE  A 

GOOD    TEST    FOR    CONVERGENCE. 

These   are   the   only   error    outputs    contained  in  the 
program    other   than   those   that    might   arise  from   library 
routines.       There   are  no    error    stops,    pauses,    dumps, 
sense   switches    or    sense   lights   used. 

Normal   output    (presuming   no   input    errors)    includes 

the  nine   sums    Sj-,    displayed   as   a   labeled   array.       Following 

this   is   one  of   the   statements: 

(SUM1) 
SUMMING  WAS    DONE    BY    SUBROUTINE    (SUM2) 

(SUM3) 

(DELTA  AND  U3  BOTH  DIFFERENT 
WHICH  ASSUMES       FROM  ZERO) 

(DELTA=0  BUT  J1+J2  AND  U3  NOT 

EQUAL  TO  ZERO) 
(Jl,  J2,  AND  DELTA  ALL  ZERO  BUT 

U3  NOT  ZERO) 

or 
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(SUM4   &   SUM6) 
SUMMING  WAS    DONE   BY    SUBROUTINES    (SUMJ+   &   SUM6) 

(SUM5   &   SUM6) 
(SUM5   &  SUM6) 

(U3    =    0    BUT    DELTA  AND    U1+U2    NOT    ZERO) 
ASSUMING    (J1,J2,    &   U3   ALL    ZERO    BUT    DELTA   NOT 

ZERO) 
(U3   AND    DELTA   ZERO   BUT    U1+U2    NOT 

ZERO) 
(Jl,    U2,    U3   AND    DELTA  ALL    ZERO) 

Following  this   is   a   statement   of   whether   or   not 

covergence  test   were   met,    and  if  not,    a   measure  of   the 

error   introduced. 
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APPENDIX   A 
PROGRAM   LISTING 


A.  IDENTIFICATION 

Title  Electric    Field   Components   by   de  Wette 

-    Schacher    Equations 
Co-op   Id.  21-NPGS-LATSUM    (F-60,63) 

Category  All   Others 

Programmer  David   Dickmann 

Date  February    1966 

B.  PURPOSE 

This   Fortram   program    computes    electric   fields   per 
unit-dipoles   for   a   lattice   of   similar    dipoles.       The 
program   accepts   lattices   of   any   symmetry   including 
triclinic   and   all  higher    symmetries,    and    computes   field 
components   in   a   cartesian   reference   frame,    using   the 
de  Wette-Schacher    equations. 

C.  USAGE 

1.  Calling   sequence:       This   is   a   main  program.       See 
input   formats. 

2.  Input   arguments: 

Al,    A2,    A3   are   the   lengths,    in  Angstroms,    of 
the   lattice   unit   cell  basis   vectors. 

XJ1,    X«J2,    XJ3   are  fractions    of   the  basis 
vectors   Al,    A2,    A3    designating   the   location  of 
the   point   at   which   the   fields   are   to   be   evaluated. 

XALPH   is   the   crystal   angle  between   the  A2   and 
A3   axes,    in  degrees   and  decimals. 

XBET    is   the   crystal   angle  between  the  Al   and 

A3   axes,    in  degrees   and   decimals. 

t 

XGAM  is   the   crystal  angle  between  the  Al   and 
A2   axes,    in   degrees   and   decimals. 
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XDEL    is    the  inclination--of   the  A3   axis   from   the 
perpendicular   to   the  A1-A2   plane,    in   degrees   and 
decimals.      If   unknown,    enter   as   XDEL=90   and 
correct   value   is   generated  by  program . 

AEX,    AEY,    AE2   are  the   absolute   error   limits 
for   the   values   of   the  x,    y,    and   z    components   of 
the   electric   fields. 

REX,    REY,    REZ   are  the   relative   error   limits 
for   the   corresponding   electric   field   components. 

N   is   a   summation   cutoff   value.       If    error   limits 
are  not    met   after  ***  4N      terms   are   summed,    the 
program    returns   the  field   values   at    that   time 
and  a   measure   of   the   maximum    error   introduced 
by   cutoff. 

3.  Space   required:      6446    cells    (F-60),    2+209    cells    (F-63) 

4.  Temporary    storage:       None 

5.  Alarms:       None 

6.  Error   returns: 

a.       Error   in   input   parameters — Control   returned 
to    monitor    with   printout: 

PROGRAM    STOPPED    BY    INPUT 

ERROR    DESIGNATED    BY    K= 

Individual   error    codes   are: 

K  Error 


1 

AU  0. 

2 

A2£  0. 

4 

A3^0. 

10 

XJ1<  0 

20 

XJ2<  0 

40 

XJT3<  0 

100 

XALPE 

200 

XBET^ 

or  XJ1>  1. 
or  XJ2>  1. 
or   XJ3>  1. 


XBET^O.    or    XBET>90. 
20 


400  XGAM^O.    or    XGAM^ISO. 

1000  XDEL<0.    or    XDEL>90. 

2000  N^O 

Multiple   errors   result   in  a    K   value   which  is 
the   sum    of   the   individual   values   above. 

b.       N    value  too    small   to    ensure   good   test   for 

convergence   of   field   sums.       Control   returned 
to    monitor    with   self-explanatory   printout. 

7.       Error    Stops:       None   except    those   associated   with 
machine   arithmetic    errors. 

8*      Input   and   output   tape    mountings:       Normal 

9.  Input   and   output   formats 

a.  Input   is    contained   on  four    cards   under   the 
following  format    control: 

READ    1,A1,A2,A3,XJ1,XJ2,XJ3,XALPH, 
XB  ET  , XGAM,  XDEL  ,  N  ,  AEX ,  R EX  , 
AEY,REY,AEZ,REZ 

1    FORMAT    (3F10.0/3F10.0/4F10.0/I2+,6E10.0) 

b.  Output:      Various 

(1)  Error   outputs   as   per   6 

(2)  S..   -  -   the  i-th   component   of   the   electric 

field  due   to   a   unit   dipole   oriented  in 
the  j-direction.       Printed   as   an  array- 
under    1PE19.10    with   labels    SUMXX,    etc. 

(3)  Hollerith   prints   indicating   subroutines 
used   and   results   of    convergence   tests. 

10.  Selective   stop   and  jump   settings:       Not   applicable 

11.  Timing: 

a.       Compile   time:      4    min   20    sec 
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b.       Run   time:      Average    1.4   sec   per   lattice, 

averaged   over  336    combinations   of   lattices, 
error   limits,    and   evaluation  points. 

12.  Accuracy:       Relative   error   of    Sn*min.    ,(REX,REY, 

RE2)*S..  J  ,J 

Absolute   errdr   of    Sns&minj   i(AEX,AEY, 
AEZ)  ,J 

13.  Equipment    configuration:       CDC    1604   with 
FORTRAN    compiler    or    CO-OP    monitor    system 
including   F-63 

14.  Library   routines   used: 

SQRTF 

SINF 

COSF 

ATANF 

ASINF 

ACOSF 

EXFF 

ABSF 

FIXFLOAT 

FLOAT  FIX 

15.  References: 

a.  D.    B.    DICKMANN,    A   Computer    Program 
for    Evaluating   Lattice   Sums,     Thesis,    U.S. 
Naval   Postgraduate   School,    Monterey,    1966 

b.  F.    W.    de  Wette   and  G.    E.    Schacher,    Phys. 
Rev.,    137,   A78    (1965) 

D.  METHOD    OF   ALGORITHM      See   references    (a)    and 
(b). 

E.  FLOW  CHART      See   reference    (a). 
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APPENDIX    B 
SUMMATION    CONVERGENCE  AND    ERRORS 

The   general   lattice   sums   to   be   evaluated   in   the 
program    are   of   the   form 


where   C.    is   a   constant,    h.    {/J, )t/J^  )    is   a  function   of /J+ 
and^6/2    which   generally   increases    with  increasing  /4*    andyd/o 
and   0[   is   a   phase  factor   of   the  form    exp  L-2.  Tt  i  JC;jUh  t/j0tyfl 
Because   of   the   phase   factors   it   is   possible   to   have  a 
distribution  of    significant   terms   over   the jJy-JjL    plane   such 
as    shown   in   figure   8. 

To    ensure   inclusion   of    significant   terms   in   such   a 
distribution   would   require    either   very   complex   convergence 
tests   or   inclusion   of    many   more   terms    then   required. 

To   avoid   this,    the   sums    were   converted   to    semi- 
infinite   sums    over   the  first    quadrant   of   the /j.-JLU    plane. 
This    led   to   the   possible   distributions    shown  in  figure  9, 
where   the   degree   of   asymmetry  about   the/f,  ~>^z   ax*s 
depended   on   the   ratio   of    the   lengths   of   the   a.,    and   a2 
basis   vectors   and   the   angle  between  them,    and   the 
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*-x/, 
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Figure  8:       Possible   distributions   of    significant   terms   of 
sum   over//.-/^   plane.      Shaded  area   contains 
all  terms   of   absolute  value   larger   than  a  given 
constant. 
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Figure   9:       Possible   distributions   of    significant   terms   of 
converted   sums 


58 


particular    sum   being    evaluated,    S-.      While   these   distributions 
still   leave   problems   in   developing   convergence   tests,    the 
problems   are  not   as    formidable   as   for   the   original    equations. 

The    solution   used   in   this    program    proceeded   in 
several   steps.       First   the   axial   terms,  /J.  =  0 ,  ^6/^=0  ,    1,    2, 
3,     .     .    and//p=0,  ^/^  =  1,    2,    3    •     •     •    "were   summed   together 
and   the   convergence   test   applied   to   the   sum   of    each   pair, 

Ml~® » Mz=1't  Ml~'x*^Z^®  *    until   the   tests    were    met   for    each 

sum    S«k,    or    N   pairs    were   summed. 

The   non-axial    terms    were   then   treated   separately. 

First   a    skew   factor,    JIG   =   2    *   MAX    (££.££)  costf+1, 

Al    A2 

was   generated.       This   number   is    such   that   the   largest   term 
of   the   sum    (in  absolute   value)    must   occur   at    some   point 
within  the  box   0£/U/,/4i^JIG   .       This    was   tested   against    N 
to    ensure   N>JJIG   so   that    the   sums    would   include   the 
largest   term. 

Summing   was    done  by   summing   all   terms  ^    =    1,    2, 
3,     .     .     •    f or  jU^-1   and   then  for^-L=2,    etc.       The    convergence 
test    were   threefold.       First   to    ensure   for   given^/V.,    that 
/J2   assumed   all   values   up   to    JIG   and   that    this   was   done 
for   all   values    of  M\    up   to    JIG.       This    ensured   picking  up 
the   most    significant   term    of   this    sum   regardless    of   its 
magnitude.       Second   to    ensure   that   for   given /A^=a.   if    soffie 

59 


term,  JA^  £.  a,    did  not    meet   the   convergence  tests   then 
summing   would  be   done  on   the  pencil /^ +=a+l, /^p~^-»    2,     .    ., 

a+1.      Last,    that   summing   would  not   exceed  the   terms 
>^l=N»/rf2=*>  Ml~U  "#2=N »    regardless   of   other   tests. 

The   flow    chart   for   this   testing  is    shown   in  figure   10. 
Actual   convergence   zones   for   a   cubic   lattice  and  a  triclinic 
lattice   are   shown   in  figures    11   and    12. 

After    completing   the  program   it   was   checked  by- 
running   sums   on    112   lattices    (distortions   of   the  basic   cubic 
lattice)    for    which   sums   are  known  from   geometric 
considerations.       The   error   analysis   for   a   relative   error 
limit    (REX,    etc)    of    10""^ ,    10"^,    and   10"8   and  an  absolute 
error   limit    (AEX,    etc)    equal  to   the  relative   limit   is   given 
in   Table   2. 
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Figure   10a      Flow    chart    of    convergence   tests    common   to 
SUM1-SUM6.       Integer   variable   I    refers   to 
values    of    summing   index  /u)  f   J      to     y(Jz. 


MAIN 
PROGRAM 


No 


Yes 


JJIG=    2.*A1*C0GAM/A2+TTI 


IendI 


fjIG=2 .  *A2*C0GAM/A1+1  .| 


1RROR 


No 


SUBROUTINE 
SUM 


■* 


Compute   terms    RTMXX,    RTMYY,    etc 
f  or  Mi  =~5-iMz  =0    and/^  =0,  MfiT^-    and 
combine 


I=I  +  l[ 


RSUMXX  =  RSUMXX+RTMXX 
etc 


Yes 


Yes 


+  REX*|  RSUMXXJt 


Form    BIGjD 
=MAX(4RTMXXf  , 
JRTMXY|  , 
etc  ) 

for   later   printout 
as    measure   of 
error   introduced 


ITMY1 
+REY* 

9 


(Continued   on  next   page) 
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Axial   terms 


Figure   10b      Continuation   of   flow    chart 
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Compute   RTMXX   =    Contribution  to    Sxx 
of   term  >/,=!  ,(U„  =  J  ,    RTMXY,    etc 
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Ye 


No 


)J=J  +  1) 
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No 

Form   LINE  1 
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of   error 
introduced 


M=l,    tag   that 
convergence   tests 
were   not    met 


No 


CONTINUE 
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+  JU, 


Figure   11: 


Actual   convergence    zones   plotted   for   unit 
cubic  lattice   using   library   subroutine 
CONTOUR. 
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X 


>A 


Figure   12: 


.Actual   convergence;    zones   plotted   for   a 
triclinic   lattice   using   library   subroutine 

CONTOUR, 
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Table   2:       Error   analysis   of   known   sums    S--  .    Ljar^eaoAftTteiw <■  the 
left-hand   numbers   are   absolute   errors,    the   right 
hand   numbers   are  relative   errors,    the  top   numbers 
are   maximum   values,    the  bottom   are   average   values, 
and   the   center   integer   is   the  number   of    sums 
evaluated. 


2a:       Input   relative   error   limit   =    10 


-4 


Input   absolute   error   limit 
All    errors   X    10~6 


=    10 


-4 


Lattice 
Type 


SUM1 


Subroutine    Used 
SUM2  SUM3 


SUM4&6        SUM5&-6 


Triclinic 

.02 

174 
.003 

.1 
.005 

0 

0 

.1 

156 

.008      . 

.1 
OOf 

0 

Monoclinic 

.03 

12 
.004 

.07 

.01 

.1 

24 
.03 

.2 

.01 

0 

.05 

8 
.008 

.09 
.01 

.04      .06 

16 
.01       .01 

Orthorhombic 

0 

-t*. 

1. 

18 
.08 

.6 
.05 

0 

0 

.7         .4 

12 
.08       .05 

Tetragonal 

0 

.2 

54 
.07 

.1 
.06 

0 

0 

.09       .1 

36 
.003    .009 

Trigonal 

0 

0 

0 

0 

0 

Cubic 

0 

.6 

18 

.07 

.1 
.01 

0 

0 

.05       .08 

12 
.005    .007 
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2b:      Input   relative   error  limit   =    10 
Input   absolute   error   limit   =    10 
All    errors   X    10~8 


-6 


-6 


Lattice 
Type 


Subroutine   Used 
SUM1  SUM2  SUM3  SUM4&6     SUM5&6 


Triclinic 

.1           1. 

174 
.02      .03 

0 

0 

1.         1. 

156 
.02        .01 

0 

Mo  no  clinic 

.2           .5 

12 
.06        .1 

1. 
.1 

24 

1. 
.1 

0 

.3            1. 

8 
.05          .2 

.2 

16 
.03 

.09 
.01 

Orthorhombic 

0 

4- 
.3 

18 

2. 
.2 

0 

0 

2. 

12 
.2 

.9 
.1 

Tetragonal 

0 

2. 

.1 

54 

1. 
.1 

0 

0 

.8 

36 
.04 

1. 
.03 

Trigonal 

0 

0 

0 

0 

0 

Cubic 

0 

6. 

18 

.9 

.07 

0 

0 

.5 

12 

.05 

.9 
.1 
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2c:      Input   relative   error   limit   =    10 
Input   absolute   error   limit   =    10 

All    errors    X    10~10 


-8 
-8 


Lattice 
Type 

SUM1 

Subroutine   Used 
SUM2          SUM3         SU-MB*»6 

SUM5&6 

Triclinic 

.2           1. 

174 
.04      .03 

0 

0 

.9           1. 
156 

.04      .02 

0 

Monoclinic 

.1           .4 

12 
.01      .04 

.9           1. 

24 
.05      .07 

0 

.1           .8 

8 
.02         .1 

.3           .08 

16 
.04      .007 

Orthorhombic 

0 

3.           1. 

18 
.2           .1 

0 

0 

2.              1. 

12 
.3             .1 

Tetragonal 

0 

3.           2. 

54 
.06      .05 

0 

0 

1.              .9 

36 
.04        .03 

Trigonal 

0 

0 

0 

0 

0 

Cubic 

0 

7.           8. 

18 
.4           -5 

0 

0 

.4             .8 

12 
.1           .07 
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APPENDIX  C 
INCOMPLETE  GAMMA  FUNCTION 

In   evaluating  the   lattice   sums   when  the   evaluation 
point   lines   in  the   x-y  plane  an  auxilliary  convergence 
function  is   used  to   smooth   and   remove  the   simple   pole   of 
the  function  being   summed.      In  the   equations   of    Schacher 
and   de  Wette  this   introduces   the  following  four   terms   into 
the   sums: 

r(s/z,7rhz) 

where  h=hCc/MyUz  )  ;  jj,  ^AJz     summing  indices;    and 

oo 

fCa.x)  S  ( ta~'e~r df  the  incomplete  gamma  function. 

In   order   to    develop   a   usable   program   based   on   Schacher 
and   de  Wette's    equation  it   was   necessary  to    obtain  a   simple 
subroutine   for    evaluating   these  functions   repeatedly  as   h 
changed  within    and^g.      Unfortunately  incomplete  gamma 
functions   are  not   widely  used   and  the  only   subroutines 
available  from   other   sources   were  those  designed  to   tabulate 
P  (a,    x)    over   a   wide   range  of   values   of   a   and  x.      Since 
P  (a,    x)    is   a   transcedential   whose    "best"   approximation  is 
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an   assymptotic    series,    tabulating   programs   are  both 
cumbersome   and   slow. 

Three    possibilities    were    considered   for   this    program. 
First,    adapt   one   of   the   existing   programs   for   limited   use. 
This    was   tried  but    even  after   abstracting   the  desired 
parts   from    the  best   of   these   programs    we   were   left    with 
a   double-precision   program    with  an   evaluation  time,    per   h 
value,    of   about   0.5    sec.       Since   we   contemplated  hundreds 
of   thousands   of    evaluations   per   lattice   we  ultimately 
discarded   this   approach. 

The   second   possibility   considered   was    simply   storing 
tabulated   values   over   the   range  h=0-100    and   interpolating 
as   necessary.       However   to    obtain  program    accuracies   of 
±10~°   for   Sj:    would  have   required   storing   4000    or    more 
values.       This,    added   to   the  large   program    storage  already 
generated   would  have   overflowed   the   internal    storage   of 
all  but   the   largest    computers   available.       Tape   or   disc 
storage   would   again  increase   run  time   too    much   so   this 
approach   was   also   discarded. 

The   final   possibility   was   adopted,    that   of   formulating 
new   approximations   to    the  incomplete   gamma   function. 
Ultimately   a   set    of    power    series   approximations    was    developed, 
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7 
using   the   Tau  Method   of   Lanczos,       to   span   the   range 

Tfn    =3-00.      Another,    long-known   series,    well    convergent 

for   small   values    of   x,    was   used  in  the  range  77/4*=0-3 

The   expansions   are 

where   the   C^ig    were   generated  using  the  Tau  Method.      The 
other   functions,  r'(3/ztx)  ,    etc   were   obtained  from  r('/z%x) 
by   using   the   recursion   relation   P  (a+/tX)  -aftOI,/)  +Xa&~*  * 
The   Tau  Method   is    well   adopted   to   exploring   several 
possibilities   of    approximating   series   on  the   computer.      In 
this   instance   we   examined   a   total   of   22+0    possibilities, 
combining   the   four   a   values,    -i,    ^,   J,   £  with  3   possible 
ranges    of   approximation,   j<X  >  I <X  }3  <X     and  in   each   case 
examining   series    of   2,    3,    h    •«•>    21    terms.       In  general, 
accuracy  increased   with   increasing  a   value,    increasing   end 
point   b,    b<x,    and   increasing   number    of   terms   in  the   series. 
We   used  a=4-   vice   f   or   I  because  it   matched  the   a   value   of 

7 
Lanczos,    C,    Applied  Analysis,    2+64-J+85,     (Prentice 

Hall,    Englewood   Cliffs,    N.J.,    1961) 
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the   expansion  for   0<x^3»      A    17    term    series    was   used 

since   it   exhausted    machine   accuracy   of   the   computer   used 

-11 
(  ±10         )  ,    and  any  higher    series    would  add   no   accuracy. 

Having   settled   on   an   approximation  it    was    checked 

against   a   double   precision   routine  known   to   be   accurate   to 

-11 
±10         .       The   graph   of   the   error   of   the   approximation  is 

shown  in  figure   13  • 
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-2x/0'« 


30 


Figure   13  s 


Error    curve   for   Incomplete  Gamma    Function 
approximation.      Abscissa   is   x,    ordinate   is 
r* V/lLtX)~r* (J/a^x)  where/**   is   approximation  and 
P  is   actual   value   of   Incomplete  Gamma   Function. 
Shaded   area   represents   an   oscillation   of    error 
of    such   a  high   frequency   that   it   could  not   be 
faithfully   represented   on   this    scale. 
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APPENDIX  D 
APPLICATION  TO  BENZENE 

During   the   preparation   of    this   program    R.M.    Hexter 
of   the  Mellon  Institute   suggested   an  application  to   a   long- 
standing  problem   in   molecular   spectroscopy.       This    work   has 
been   completed   and  is    detailed  below   as   an   example   of   applying 
the   program   to   an  actual  problem. 

In   calculating  the   electronic    spectrum   of   benzene   it 
is   necessary   to    compute   the  transition   moments    of   the   free 
molecules.       These   moments   depend   on  the   dipole-dipole   sums 
of   the   crystal.      Various    workers    computing   these   sums   by 

various    methods   have   obtained   varying  values   for   the   sums. 

8 
Craig   and  Walsh      have   set   out   the   crystalline   structure 

assumed  and  tabulated   the   values    obtained  by  two    sunmation 

procedures.       This   program    was   used   to    check   these   values 

with    results   tabulated  below. 

The  benzene   crystal    model   consists    of   an   orthorhombic 

unit    cell   with   four    molecular   orientations   at   the   corners 

and  faces.       If     we   superimpose   cartesian   reference   frames 

I 
on   each    molecule   we   can  identify  our    molecules   as   follows: 


8D.    Craig   and   J.    Walsh,    J.    Chem    Soc    1958,    1613 


Molecular       Unit   Cell  Crystal  Direction   Cosines 

Type  Location  Axis  x  y 

I  0,0,0  a 

b 
c 


11  i>°>i 


a 
b 
c 


.6365 

.2000 
-.7450 

-.2992 

.9541 
0 

.7109 
.2229 
.6671 

-.6365 
-.2000 
-.7450 

.2992 

-.9541 
0 

-.7109 

-.2229 

.6671 

-.6365 
.2000 
.7450 

.2992 
.9541 

0 

-.7109 

-.2229 

.6671 

.6365 

-.2000 

.7450 

-.2992 
-.9541 

0 

.7109 
-.2229 
-.6671 

III  0,4,4  a 

b 

c 

IV  4,4,0  a 

-*  b 

c 

Where    (a,b,c)    forms   the   cartesian  frame   of   the  unit   cell. 
The  unit   cell   parameters   are    (at   -3°    C) 

a-L    =    7.460   A  cr  =/3"ar  =  90° 

a2    =   9.666  A 

a3    =   7.034  A 

The   problem    then  is   to    compute   the  interaction   energy 
of   a  unit   dipole   in  the   i-   direction  in   molecule  I    due  to   a 
lattice     of   unit   dipoles   in  the  j-direction  in   molecules   of 
type    >•  ,    where   i,    j    =   x,    y,    z;    V  =   I,    II,    III,    IV. 

In   comparing  figures   obtained  by   this   program    with 
those   tabulated  by   Craig   and  Walsh   two   difficulties   arise. 
First   dipole-dipole   energies   from   this   program    work   out  in 
recipricol  angstroms   per   coulomb   squared   while   Craig   and 
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Walsh   tabulate   in   recipricol   centimeters    (wave   number)    per 
electronic    charge    squared.       This    requires    multiplying   the 
results   of    this   program   by    1.161    X    10      to   be    comparable. 

Second   Craig   and  Walsh   assume   an  infinite   spherical 
crystal    while   this   program   assumes   an  infinite   slab-shaped 
crystal.       Thus    depolarization   factors    must   be   applied.       The 
two   systems   are: 

Sperical  Slab-shaped 

Thus,    to   be   comparable    we   must   add   -    {k-Tf  /vo\)    = 
0.0165170    to    23     terms   generated  by   this   program   and 
subtract   -    {l+V  /vol)    =   0.0082585   to   xx  and   xy   terms 
generated  by   this   program.       (Where   vol   =   volume   of   unit 
cell  in  A) 

GENERAL    PROCEDURE 

With  these   considerations   in   mind   the   computation 
becomes    straight-forward  if   tedious.       First    we   will   compute 
electric   field   components   at    molecule   I    along   the   crystal 
axes   a,    b,    c   due   to   an   assumed   set   of   unit   dipoles   along 
the   a,    b,    and   c   axes   at    molecular   positions   I,    II,    III,    and 
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* 
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IV.       This   is   done  by   calculating   the  field  at   the  appropriate 

position  of    molecule   I    (given  by  XJI,    XJ2,    XJ3)    within 

each   of   the  I,    II,    III,    and  IV    sublattice  unit   cells.      Then 

we   correct   the  depolarization  factors   and  apply  the   conversion 

factor.      Then,    using  the  direction   cosines   of   the  type  I 

molecule  for    each   molecular   axis   x,    y,    and   z    (with  respect 

to    the  a,    b,    c   axes)    we   can  get   the  x,    y,    and   z   components 

of   field  at   I    due  to   a,    b,    and   c   dipoles   at   all   molecular 

positions.       Then  by   multiplying  by   the   direction   cosines   for 

each   molecule   respectively   we   get    energies    of   the  x,    y,    and 

z    dipoles    of    molecule   I    due   to   x,    y,    and   z    dipoles   at    each 

molecular   position.       These   we   compare   with   similarly 

tabulated  results    of    Craig  and  Walsh. 

DETAILS    OF    COMPUTATION 

To   get   the  a,    b,    and   c   field   components   at    molecule 

I    we   enter   the   program    with: 

Al   =   7.460  XALPH   =   90. 

A2    =   9.666  XBET    =   90. 

A3    =    7.034  XGAM   =    90. 

XDEL    =      0. 

AEX=REX=AEY=REY=AE2=REZ=10"7 

N=20 

and  for   various    molecules: 
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II 


III 


IV 


XJ1 

0 

* 

0 

ft 

XJ2 

0 

0 

ft 

ft 

XJ3 

0 

* 

ft 

0 

Since   the   crystal   is   orthorhombic    we   get    S..=0,    i^j. 
The   three   nonzero    outputs    SUMXX,    SUMYY ,    and   SUM2Z 
here   correspond   to   the   a,    b,    and   c   field   components.       The 
computed   Sj,-   values,    depolarization   correction,    conversion 
factor,    and   converted   values   are   listed  below: 


Fields   at  i-field 

I    due   to  component 

molecules  due  to  j-  Computed 

of   type  dipoles  Values 


Corrected 
Depolarization        Values 


II 


III 


IV 


aa 
bb 
cc 

.0104233 
.0011720 

.0115953 

-.0082585 

-.0082585 

.0165170 

.0021648 
.0070865 
.0049217 

aa 
bb 
cc 

.0240220 
-.0187985 
-.0052235 

-.0082585 

-.0082585 

.0165170 

.0157635 

-.0270570 

.0112935 

aa 
bb 
cc 

-.0055470 

.0231906 

-.0176436 

-.0082585 

-.0082585 

.0165170 

-.OI3-8O  5  5 

.0149321 

-.0011266 

aa 
bb 
cc 

.0080050 

.0209184 

-.0289234 

-.0082585 

-.0082585 

.0165170 

-.0002535 

.0126599 

-.0124064 
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Molecule 

Axes 

Corrected 
Value 

Conversion 
Factor 

Converted 
Values 

I 

aa 
bb 
cc 

.0021648 

-.0070865 

.0049217 

1.161X105 

251.333 

-822.742 

571.409 

II 

aa 
bb 
cc 

.0157635 

-.0270570 

.0112935 

1830. 142 

-3141.317 

1311.175 

III 

aa 
bb 
cc 

-.0138055 

.0149321 

-.0011266 

-1602.818 
1733. 616 
-130.798 

IV 

aa 
bb 
cc 

-.0002535 

.0126599 

-.0124064 

1.161X105 

-29.431 

1469. 814 

-1440.383 

Now   we  need  the   x,    y,    and   z   components   of   these 
fields   at    molecule  I;    I-I    (x,  m),    I-I    (y,  m),    I-I    (z,  m), 
etc    (m=a,b,c).       Here   we  need   the   molecule   I    direction 
cosines   for   x,    y,    or    z   with  respect   to   a,    b,    and   c*      Two 
examples    suffice: 

x     -   field    components 


Molecular 
Dir.    Cos      Component 


Fields 
due  to 
at 

at   I 
dipoles 

Unit    Call 
Component 

Value 

I 

aa 
bb 
cc 

251.333 

-822.742 

571.409 

II 

1 

1 
t 

1 

1 

.6365       159.9734 

.2000    -164.5484 
-.7450    -425.6997 
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y    -   field   components 

Fields   at   I 

due  to   dipoles       Unit    Call  Molecular 

at  Component        Value  Dir   Cos       Component 

i 
i 
i 

II  aa  1830.12+2    -.2992        -547-5784 

bb  -3141.317       .9541        -2997.1305 

cc  1311.175  0  0 

III 

1 

1 

From    these   we   can   make  up   the   energies   by   multiplying 
the    molecular    components   by   the  direction   cosines    of   the 
source   dipole   at   the   molecular   position  of   interest. 

Since     we   have   tabulated  the   molecular    components 
I-I    (x,    m)    and   I-II    (y,  m),    if   we   want   the  energy  of   an 
x-dipole  in   molecule  I    due  to   a   lattice  of   y-dipoles   at 
position   I,    I-I    (x,y),    we   need   only   multiply   the    molecular 
components   I-I    (x,    m)    by   the  y   direction   cosines   of   type 
I    and   sum.      We   show   the   computation  of   I-I    (x,y)    and 
I-II    (y,z)    below: 
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Field                                                         Direction  Energy 

Component        Axis       Value        Cosine  Contribution 

Iy   dipoles 

I-I    (x,       )           aa      159-9734      -.2992  -47.8640 

bb  -164.5484        -9541  -156.9956 

cc  -425.6997            0  0 

I-I    (x,y)    energy -204.8596 

IIz    dipoles 
I-II    (y,      )        aa  -547-5784     -.7109  389.2734 

bb-2997.1305     -.2229  668.0603 

cc  0  .6671  0 

I-II    (y,z)    energy  1057.3337 

These   then  are   the  numbers   to   be   compared   with 

those   tabulated  by   Craig   and  Walsh.       Their   table   contains 

Q 

numbers   obtained  both   from   the   Ewald-Kornf eld      method 
applied   to   an  infinite      sphere   and  from   a   direct    molecule- 
molecule   summation  over   a   sphere   of    radius   20   A.       The 
comparison  is   shown  below.      Results    (a)    are  from   this 
program,    (b)    from   the    Ewald-Kornf eld   method,    and    (c) 
from    direct    molecule-molecule   summation. 


9 
Kittel,    Introduction  to    Solid   State   Physics,    2nd    Ed, 

347,    (vJohn  Wiley  and   Sons,    Inc.,    New    York,    1954) 
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Dipole  Molecule  Sites 

Axes           Results        I,    I                I,    II  I,    III  I,    IV 

x,    x               (a)  -386.06        -111.94  -791.29  -728.72 

(b)  -386               -112  -791  -729 

(c)  -415                -92  -769  -771 

y,  y      (a)  726.45  -2695-73  -1721.61  1340.62 

(b)  727  -2679  -1723  1341 

(c)  790  -2629  -1713  1405 

z,  z      (a)  -340.47    185-49  -954-10  -553-20 

(b)  -341       185  -955  -553 

(c)  -375       220  -950  -587 

x,  y      (a)  204-86  -947-96  -25-57  274-87 

(b)  205  -947  -25  275 

(c)  233  -934  -8  289 

x,  z      (a)  206.94   1339-72  -737-54  794-70 

(b)  207      1340  -738  795 
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